!-------------------------------------------------------------------------------
! FileName:     gol_mpi_nbym_parallel_io.F90
! 
! Compile :     $ mpif90 -g -o gol_mpi_nbym_parallel_io.x -Wall gol_mpi_nbym_parallel_io.F90
!
! Run     :     $ mpirun -np 4 -hostfile hosts gol_mpi_nbym_parallel_io.x
!
! Be careful:   1. GridX, GridY must be same size. (ex) 10x10, 20x20, 100x100
!               2. To do preprocessor, .F90 is required.
!               3. # of processor is 4, 9, 16, 25....f
!
! Tab size:     8
! Version :     1.03
!
! Updated by    2011. 11. 03
!-------------------------------------------------------------------------------


PROGRAM gol_mpi_nbym_parallel_io
IMPLICIT NONE
INCLUDE 'mpif.h'
        INTEGER :: ROOT = 0
        INTEGER nRank, nProcs, iErr, status(MPI_STATUS_SIZE)
        INTEGER, PARAMETER ::   nGridXSize=20, nGridYSize=20, nSteps=100, nSleep=1
        INTEGER nSqrtNprocs, nMyWidthX, nMyWidthY, nLeftOver, nQuotient, nRemainder
        INTEGER nMyStartX, nMySTartY
        INTEGER nTmpMyStartX, nTmpMyStartY, nTmpMyWidthX, nTmpMyWidthY
        INTEGER i, j, k, l
        INTEGER nLiveCell
        INTEGER mpiFile
        
        INTEGER, DIMENSION(:, :), ALLOCATABLE                   :: gol_buf
        INTEGER, DIMENSION(:, :), ALLOCATABLE                   :: tmp_gol_buf
        
        INTEGER column_new_type
        
        ! MPI Start     
        CALL MPI_INIT(iErr)
        CALL MPI_COMM_SIZE(MPI_COMM_WORLD, nProcs, iErr)
        CALL MPI_COMM_RANK(MPI_COMM_WORLD, nRank, iErr)
        
        ! Let each processor decide what column(s) to run       
        IF (nProcs == 2) THEN
                nSqrtNprocs = 2
        ELSE
                nSqrtNprocs = SQRT(REAL(nProcs))
        ENDIF
        
        nMyWidthX = nGridXSize / nSqrtNprocs
        nMyWidthY = nMyWidthX
        nLeftOver = MOD(nGridXSize, nSqrtNprocs)
        
        nQuotient = nRank / nSqrtNprocs
        nRemainder = MOD(nRank, nSqrtNprocs)
        
        nMyStartX = nQuotient * nMyWidthY
        nMyStartY = nRemainder * nMyWidthY
        
        ! Right edge
        IF (MOD(nRank+1, nSqrtNprocs) == 0) THEN
                nMyWidthY = nMyWidthY + nLeftOver
        ENDIF
        
        ! Bottom edge
        IF (nProcs == 2) THEN
                nMyWidthX = nGridYSize
        ELSEIF ((nRank/nSqrtNprocs) == (nSqrtNprocs-1)) THEN
                nMyWidthX = nMyWidthX + nLeftOver
        ENDIF

        ALLOCATE(gol_buf(0:nGridXSize, 0:nGridYSize))
        
        ! set to 0 
        DO i=0, nGridXSize-1
                DO j=0, nGridYSize-1
                        gol_buf(i,j) = 0
                END DO
        END DO
        
        if (nRank == ROOT) THEN
                WRITE (*, *) 'nProcs = ', nProcs
                WRITE (*, *) 'nRank = ', nRank
                WRITE (*, *) 'nSqrtNprocs = ', nSqrtNprocs
                WRITE (*, *) 'nMyWidthX = ', nMyWidthX
                WRITE (*, *) 'nMyWidthY = ', nMyWidthY
                WRITE (*, *) 'nLeftOver = ', nLeftOver
        
                CALL getRandomValue(gol_buf, nGridXSize, nGridYSize)
                
                !DO i=0, nGridXSize-1
                !       DO j=0, nGridYSize-1
                !       PRINT *, i, j, gol_buf(i,j)
                !END DO
                
                CALL showGridState(gol_buf, nGridXSize, nGridYSize)
                WRITE(*, *)     
        ENDIF
        
        ! Allocate memory buffer for holding own GOL values
        IF (nRank .NE. ROOT) THEN
                ALLOCATE(tmp_gol_buf(0:nMyWidthX, 0:nMyWidthY))
        END IF

        ! ----------------------------------------------------------------------
        ! Main Loop : time steps
        ! ----------------------------------------------------------------------
        DO i=1, nSteps
                CALL MPI_BARRIER(MPI_COMM_WORLD, iErr)

                ! Let all processor have same data              
                CALL MPI_TYPE_CONTIGUOUS(nGridXSize, MPI_INTEGER, column_new_type, iErr)
                CALL MPI_TYPE_COMMIT(column_new_type, iErr)
                
                DO j=0, nGridYSize-1
                        CALL MPI_BCAST(gol_buf(0, j), 1, column_new_type, ROOT, MPI_COMM_WORLD, iErr)
                        !CALL MPI_BCAST(gol_buf(0, j), nGridXSize, MPI_INTEGER, ROOT, MPI_COMM_WORLD, iErr)
                END DO

                !CALL MPI_TYPE_FREE(column_new_type)

                ! Do work at the my own grid cell
                CALL calEvolution(gol_buf, nGridXSize, nGridYSize, nMyStartX, nMyStartY, nMyWidthX, nMyWidthY)

                ! Send my own do-work data to ROOT(0)
                IF (nRank .NE. ROOT) THEN
                        DO j=nMyStartX, nMyStartX+nMyWidthX-1
                                DO k=nMyStartY, nMyStartY+nMyWidthY-1
                                        tmp_gol_buf(j-nMyStartX, k-nMyStartY) = gol_buf(j, k)
                                END DO
                        END DO

                        CALL MPI_SEND(nMyStartX, 1, MPI_INTEGER, ROOT, 20, MPI_COMM_WORLD, iErr)
                        CALL MPI_SEND(nMyStartY, 1, MPI_INTEGER, ROOT, 21, MPI_COMM_WORLD, iErr)
                        CALL MPI_SEND(nMyWidthX, 1, MPI_INTEGER, ROOT, 22, MPI_COMM_WORLD, iErr)
                        CALL MPI_SEND(nMyWidthY, 1, MPI_INTEGER, ROOT, 23, MPI_COMM_WORLD, iErr)

                        CALL MPI_TYPE_CONTIGUOUS(nMyWidthX, MPI_INTEGER, column_new_type, iErr)
                        CALL MPI_TYPE_COMMIT(column_new_type, iErr)
                        
                        DO j=nMyStartY, nMyStartY+nMyWidthY-1
                                CALL MPI_SEND(tmp_gol_buf(0, j-nMyStartY), 1, column_new_type, ROOT, 30, MPI_COMM_WORLD, iErr)
                        END DO

                        !CALL MPI_TYPE_FREE(column_new_type)

                ELSE IF (nRank == ROOT) THEN
                        DO j=1, nProcs-1
                                CALL MPI_RECV(nTmpMyStartX, 1, MPI_INTEGER, j, 20, MPI_COMM_WORLD, status, iErr)
                                CALL MPI_RECV(nTmpMyStartY, 1, MPI_INTEGER, j, 21, MPI_COMM_WORLD, status, iErr)
                                CALL MPI_RECV(nTmpMyWidthX, 1, MPI_INTEGER, j, 22, MPI_COMM_WORLD, status, iErr)
                                CALL MPI_RECV(nTmpMyWidthY, 1, MPI_INTEGER, j, 23, MPI_COMM_WORLD, status, iErr)
                                
                                ! WRITE(*, *) 'nRank=', nRank, 'nTmpMyStartX=', nTmpMyStartX, 'nTmpMyStartY=', nTmpMyStartY, 'nTmpMyWidthX=',&
                                !       nTmpMyWidthX, 'nTmpMyWidthY=', nTmpMyWidthY
                                                        
                                ALLOCATE(tmp_gol_buf(0:nTmpMyWidthX, 0:nTmpMyWidthY))

                                CALL MPI_TYPE_CONTIGUOUS(nTmpMyWidthX, MPI_INTEGER, column_new_type, iErr)
                                CALL MPI_TYPE_COMMIT(column_new_type, iErr)
                                
                                DO k=0, nTmpMyWidthY-1
                                        CALL MPI_RECV(tmp_gol_buf(0, k), 1, column_new_type, j, 30, MPI_COMM_WORLD, status, iErr)
                                END DO
                                
                                !CALL MPI_TYPE_FREE(column_new_type)

                                DO k=nTmpMyStartX, nTmpMyStartX+nTmpMyWidthX-1
                                        DO l=nTmpMyStartY, nTmpMyStartY+nTmpMyWidthY-1
                                                gol_buf(k, l) = tmp_gol_buf(k-nTmpMyStartX, l-nTmpMyStartY)
                                        END DO
                                END DO
                        
                                DEALLOCATE(tmp_gol_buf)
                        END DO
                END IF

                CALL MPI_BARRIER(MPI_COMM_WORLD, iErr)
                
                IF (nRank == ROOT) THEN
                        WRITE(*, *)
                        WRITE(*, '(A, I6, A, I6)') 'msg> nStep = ', i, '/', nSteps
                        WRITE(*, *)
                        CALL showGridState(gol_buf, nGridXSize, nGridYSize)
        
                        nLiveCell = 0
        
                        DO j=0, nGridXSize-1
                                DO k=0, nGridYSize-1
                                        nLiveCell = nLiveCell + gol_buf(j, k)
                                END DO
                        END DO
                        
                        WRITE(*, *)
                        WRITE(*, '(A, I4)') 'msg> nLiveCell = ', nLiveCell
        
                        WRITE(*, '(A, I2, A)') 'msg> Wait a ', nSleep, ' sec'
                        CALL Sleep(nSleep)
                END IF
        END DO
        ! ----------------------------------------------------------------------
        ! End of Main Loop : time steps
        ! ----------------------------------------------------------------------
        
        IF (nRank .NE. ROOT) THEN
                DEALLOCATE(tmp_gol_buf)
        END IF
        
        ! MPI Parallel Io
        IF (nRank == ROOT) THEN
                CALL MPI_FILE_OPEN(MPI_COMM_SELF, 'mpi_parallel_io_out_fortran.txt', MPI_MODE_WRONLY+MPI_MODE_CREATE, &
                                        MPI_INFO_NULL, mpiFile, iErr)
                CALL MPI_FILE_SET_VIEW(mpiFile, 0, MPI_INTEGER, MPI_INTEGER, 'native', MPI_INFO_NULL, iErr)     

                DO i=0, nGridXSize-1
                        DO j=0, nGridYSize-1
                                CALL MPI_FILE_WRITE(mpiFile, gol_buf(i, j), 1, MPI_INTEGER, MPI_STATUS_IGNORE, iErr)
                        END DO
                END DO

                CALL MPI_FILE_CLOSE(mpiFile, iErr)
        END IF


        DEALLOCATE(gol_buf)

        ! MPI End
        CALL MPI_FINALIZE(iErr)

!------------------------------------------------------------
CONTAINS
SUBROUTINE getRandomValue(gol_buf, nGridXSize, nGridySize)
IMPLICIT NONE
INTEGER nGridXSize, nGridYSize
INTEGER, DIMENSION(0:nGridXSize, 0:nGridYSize)          :: gol_buf

        INTEGER ii, jj
        REAL    arand
        
        CALL INIT_RANDOM_SEED()
        
        DO ii=0, nGridXSize-1
                DO jj=0, nGridYSize-1
                        CALL RANDOM_NUMBER(arand)
                        gol_buf(ii, jj) = NINT(arand)           ! Get random # (0 or 1) 
                END DO  
        END DO
END SUBROUTINE

SUBROUTINE showGridState(gol_buf, nGridXSize, nGridYSize)
IMPLICIT NONE
INTEGER nGridXSize, nGridYSize
INTEGER, DIMENSION(0:nGridXSize, 0:nGridYSize)          :: gol_buf
        INTEGER ii, jj
        
        WRITE(*, 400, advance='no')     '   '
        
        DO ii=0, nGridYSize-1
                WRITE(*, 100, advance='no')     ii
        END DO
        
        WRITE(*, *)                             ! new line
        WRITE(*, 400, advance='no')     '   '
        
        Do ii=0, nGridYSize-1
                WRITE(*, 200, advance='no')     '--'
        END DO
        
        WRITE(*, *)                             ! new line
                
100     FORMAT(I2, X)
200     FORMAT(A2, X)

        DO ii=0, nGridXSize-1
                WRITE(*, 300, advance='no')     ii, '|'
                
                DO jj=0, nGridYSize-1
                        IF (gol_buf(ii, jj) == 0) THEN
                                WRITE(*, 400, advance='no')     ' - '
                        ELSE IF (gol_buf(ii, jj) == 1) THEN
                                WRITE(*, 400, advance='no')     ' # '
                        ELSE
                                WRITE(*, 400, advance='no') ' $ '
                        END IF
                END DO
                WRITE(*, *)                             ! new line
        END DO 

300     FORMAT(I2, A)
400     FORMAT(A)
END SUBROUTINE

FUNCTION getNeighborCount(ii, jj, tmp_buf, nGridXSize, nGridYSize)
INTEGER ii, jj, nGridXSize, nGridYSize
INTEGER, DIMENSION(0:nGridXSize, 0:nGridYSize)          :: tmp_buf
INTEGER getNeighborCount
        INTEGER a, b, c, d, e, f, x
        INTEGER count

        a = jj - 1              
        b = jj                  
        c = jj + 1              
        d = ii - 1              
        e = ii                  
        f = ii + 1              
        count = 0
        
        DO x=0, 7
                SELECT CASE(x)
                        CASE (0)
                                IF (a >=0 .AND. d >= 0)                 count = count + tmp_buf(a, d)
                        CASE (1)
                                IF (a >= 0)                             count = count + tmp_buf(a, e)
                        CASE (2)
                                IF (a >= 0 .AND. f < nGridYSize)        count = count + tmp_buf(a, f)
                        CASE (3)
                                IF (f < nGridYSize)                     count = count + tmp_buf(b, f)
                        CASE (4)
                                IF (c < nGridXSize .AND. f < nGridYSize)count = count + tmp_buf(c, f)
                        CASE (5)
                                IF (c < nGridXSize)                     count = count + tmp_buf(c, e)
                        CASE (6)
                                IF (c < nGridXSize .AND. d >= 0)        count = count + tmp_buf(c, d)
                        CASE (7)
                                IF (d >= 0)                             count = count + tmp_buf(b, d)
                END SELECT
        END DO  

        getNeighborCount = count
END FUNCTION

SUBROUTINE calEvolution(gol_buf, nGridXSize, nGridYSize, nMyStartX, nMyStartY, nMyWidthX, nMyWidthY)
INTEGER nGridXSize, nGridYSize, nMyStartX, nMyStartY, nMyWidthX, nMyWidthY
INTEGER, DIMENSION(0:nGridXSize, 0:nGridYSize)                  :: gol_buf
        INTEGER ii, jj
        INTEGER, DIMENSION(0:nGridXSize, 0:nGridYSize)          :: tmp_buf
        INTEGER nNeighbors
        
        ! Copy from gol-buf to tmp_buf
        ! DO ii=0, nGridXSize-1
        !       DO jj=0, nGridYSize-1
        !               tmp_buf(ii, jj) = gol_buf(ii, jj)
        !       END DO
        ! END DO

        ! Copy from gol-buf to tmp_buf
        tmp_buf(0:nGridXSize, 0:nGridYSize) = gol_buf(0:nGridXSize, 0:nGridYSize)
        
        DO ii=nMyStartX, nMyStartX+nMyWidthX-1
                DO jj=nMyStartY, nMyStartY+nMyWidthY-1
                        nNeighbors = getNeighborCount(jj, ii, tmp_buf, nGridXSize, nGridYSize)
                        
                        IF (nNeighbors == 2) THEN               ! no change
                                ! gol_buf(ii, jj) = gol_buf(ii, jj)
                        ELSE IF (nNeighbors == 3) THEN          ! alive
                                gol_buf(ii, jj) = 1
                        ELSE                                    ! dead
                                gol_buf(ii, jj) = 0
                        END IF
                END DO
        END DO
END SUBROUTINE

SUBROUTINE INIT_RANDOM_SEED()
            INTEGER :: i, n, clock
            INTEGER, DIMENSION(:), ALLOCATABLE :: seed
          
            CALL RANDOM_SEED(size = n)
            ALLOCATE(seed(n))
          
            CALL SYSTEM_CLOCK(COUNT=clock)
          
            seed = clock + 37 * (/ (i - 1, i = 1, n) /)
            CALL RANDOM_SEED(PUT = seed)
          
            DEALLOCATE(seed)
END SUBROUTINE

!-------------------------------------------------------------------------------
END
